version 12.1
set more off
program drop _all
  quietly log
  local logon = r(status)
  if "`logon'" == "on" { 
	log close 
	}
log using bds2020psrm.log, text replace

/*		********************************************************************	*/
/*     	File Name:		bds2020psrm.do											*/
/*     	Date:   		February 28, 2020										*/
/*      Author: 		Frederick J. Boehmke									*/
/*      Purpose:		Replicate analysis of filibuster and crisis duration	*/
/*						duration reported in Tables 1 and 2 and Figure 2. 	 	*/
/*						in Boehmke, Dion, and Shipan (2020, PSRM). 				*/
/*      Input File:		bds2020psrm-filibuster.dta								*/
/*						bds2020psrm-crises										*/
/*      Output File:	bds2020psrm.log											*/
/*						bds2020psrm-filib-alltstar.txt							*/
/*						bds2020psrm-filib-loglik.dta							*/
/*						bds2020psrm-TXMY.ster									*/
/*						bds2020psrm-table01.txt									*/
/*						bds2020psrm-table02.txt									*/
/*						bds2020psrm-figure02.gph								*/
/*						bds2020psrm-figure02.gph								*/
/*		Requires:		warofatt.ado, 											*/
/*						estout.ado												*/
/*		********************************************************************	*/

	/**********************************************/
	/* Uncomment to install needed Stata add-ons. */
	/**********************************************/

	
*ssc install estout
*net install warofatt, from(https://myweb.uiowa.edu/fboehmke/stata/)
	

	/********************************************************************************/
	/* This provides a Weibull MLE for proper likelihood ratio tests with warofatt. */
	/********************************************************************************/

	
program define wblreg

	args lnf theta1 theta2

	tempvar lambda1 durdep

	quietly gen double `lambda1'	= exp(`theta1')
	quietly gen double `durdep' 	= exp(`theta2')
	
	quietly replace `lnf' = `theta1' + ln(`durdep') + (`durdep'-1)*ln(`lambda1'*$ML_y1) ///
				- (`lambda1'*$ML_y1)^`durdep'

end

	
	/**************************************/
	/**************************************/
	/* Do the analysis of the filibuster. */
	/**************************************/
	/**************************************/
	
		/* Read in the replication data from Dion, Boehmke, MacMillan and Shipan. */

use bds2020psrm-filibuster, clear
	
		/* Label the variables we use to make the table look nicer. */
	
	label variable days 			"Filibuster Length"
	label variable yesbinder 		"High Importance Policy"
	label variable unclearbinder 	"Unclear Importance Policy"
	label variable post1975			"1875 and After"

	
	/***************************************/
	/* Run the models reported in Table 1. */
	/***************************************/

	
		/* Start with the basic Weibull models. Use our likelihood to get */
		/* comparable log-likelihood values for ests between Weibull and warofatt.  */
	
		/* Model 1. */
	
ml model lf wblreg (days = ) /ln_p, robust

	ml maximize

		/* Calculate the adjusted AIC to add to the table. */

	estadd scalar aic_adj = -2*e(ll) + 2*e(k)

	estimates store T1M1
	estimates save bds2020psrm-T1M1, replace

		/* Model 2. */
	
ml model lf wblreg (days = yesbinder unclearbinder post1975) /ln_p, robust

	ml maximize

	estadd scalar aic_adj = -2*e(ll) + 2*e(k)

	estimates store T1M2

		/* Model 3: Specify the war of attrition model without covariates. */
		/* This plots the hazard and saves the log-likelihood results for */
		/* each value of t* so we can plot them later. */
	
warofatt days, robust startvalues iterate(50) difficult ///
	saveloglik(bds2020psrm-filib-loglik, replace)
	
	estadd scalar aic_adj = -2*e(ll) + 2*(e(k) + 1)

	estimates store T1M3
	estimates save bds2020psrm-T1M3, replace
	
		/* Model 4: Now add covariates in two equations. */
	
warofatt days yesbinder unclearbinder post1975, probss(yesbinder unclearbinder post1975) /// 
	robust startvalues iterate(50) difficult

	estimates store T1M4
	estimates save bds2020psrm-T1M4, replace
  
	estadd scalar aic_adj = -2*e(ll) + 2*(e(k) + 1)
	
		/* Test equality between policy importance coefficients by equation. */

	lincom [Weak]_b[yesbinder] - [Weak]_b[unclearbinder]
	lincom [logit_pi]_b[yesbinder] - [logit_pi]_b[unclearbinder]
	
	
  	/**********************************************************/
	/* Do chi2 tests for model fit.	Note that the critical    */
	/* value needs to be adjusted via the chi-bar correction. */
  	/**********************************************************/


lrtest T1M3 T1M1, force
lrtest T1M4 T1M2, force

   
	/************************************************************/
	/* This runs the model separately for each value of t* in   */
	/* order to save those results to compare parameter values.	*/
	/************************************************************/

		/* Fractional values of tstar work, but add nothing since */
		/* classification is based on integer realizations in the data. */
  
levelsof days, local(levdays)

foreach tstar of local levdays {

  capture warofatt days, tstarvalues(`tstar') startvalues iterate(50) difficult

  if e(converged) == 1 {

	display in green `tstar'
	
	quietly estadd scalar aic_adj = -2*e(ll) + 2*(e(k) + 1)
	
	estimates store tstar_`tstar'_all
	  
	}

  else {
  
	display in red `tstar'

	}
	
  }

    
  	/*****************************************************************/
	/* Create table 1 and a table comparing across all values of t*. */
  	/*****************************************************************/


estout T1M1 T1M2 T1M3 T1M4 ///
	using bds2020psrm-table01.txt, ///
	cells(b(star fmt(3)) se(par)) ///
	starlevels(* 0.10 ** 0.05) ///
	equations(1:1:1:1 , 2:2:4:4) ///
	order(#1: #2: StrongStrong: ln_p_SS: logit_pi:) ///
	eqlabels("Weak Hazard" "Weak Dur. Dep." "SS Hazard" "SS Dur. Dep." "Probability Cured") ///
	stats(N ll tstarmax Fwtstar aic_adj, fmt(0 2 0 4 2 2)) ///
	label modelwidth(11) varwidth(30) ///
	varlabels(_cons constant) ///
	style(tab) replace
	
estout tstar_*_all ///
	using bds2020psrm-filib-alltstar.txt, ///
	cells(b(star fmt(2)) se(par) ) starlevels(* 0.10 ** 0.05) ///
	stats(N ll tstarmax Fwtstar aic aic_adj, fmt(0 2 0 4 2 2)) ///
	label modelwidth(6) varwidth(30) ///
	eqlabels(Weibull "Duration Dependence" Exponential "Probability Cured") ///
	varlabels(_cons constant) ///
	style(tab) replace

	
  	/**********************************************/
	/* Figure 2: Graph the final log-likelihoods. */
  	/**********************************************/


use bds2020psrm-filib-loglik, clear

		/* Identify the maximum value of the reported */
		/* log-likelihoods across all candidates values of t*. */

  quietly sum ll_value if ll_convg == 1 & ll_iters < ll_maxiters

		/* Now find the smallest value of t* that produces the maximum LL. */
  
	quietly sum ll_tstar if ll_value == r(max)
	local tstar_max = r(min)
	levelsof ll_tstar if ll_conv!=1 & ll_iters==ll_maxiters & ll_tstar<50
	
  twoway connected ll_value ll_tstar if ll_conv==1 & ll_iters<ll_maxiters & ll_tstar<50, scheme(s1mono) ///
	lcolor(gs0) mcolor(gs0) ///
	xtitle(t*) ///
	xlabel(0 `tstar_max' 20 30 40 50) ///
	xtick(0(10)50) ///
	xtick(`r(levels)', tposition(inside) tlwidth(medthin) tlength(*1.5)) ///
	ytitle(Final Log-likelihood) ///
	ylabel(#6, grid) ///
	xline(`tstar_max', lpattern(dash) lwidth(thin)) ///
	saving(bds2020psrm-figure02, replace)


	/**************************************/
	/**************************************/
	/* Do the analysis of crisis length.  */
	/**************************************/
	/**************************************/

		/* Read in the replication data from DeRouen and Goldfinch. */

import excel using bds2020psrm-crises.xls, clear firstrow case(lower)

	generat logrelcap = ln(relcap)

		/* label the variables we use to make the table look nicer. */
	
	label variable sevvio		"Violence"
	label variable trgterra  	"Crisis Duration"
	label variable logrelcap  	"Log Relative Cap."
	label variable contig  		"Contiguity"
	label variable ethnic  		"Ethnic"
	label variable terr  		"Territory"
	label variable dyaddem3  	"Joint Democracy"
	label variable socun  		"Social Unrest"
	label variable rival 		"Rivals"


	/************************************/
	/* Run the various specifications.	*/
	/************************************/

		/* Start with the basic Weibull models. Use our likelihood to get */
		/* comparable log-likelihood values for ests between Weibull and warofatt.  */
	
		/* Model 1. */
	
ml model lf wblreg (trgterra = ) /ln_p, robust

	ml maximize

	estadd scalar aic_adj = -2*e(ll) + 2*e(k)

	estimates store T2M1
	estimates save bds2020psrm-T2M1, replace
    
ml model lf wblreg (trgterra =  sevvio logrelcap dyaddem3 contig rival ethnic terr socun) /ln_p, robust

	ml maximize

	estadd scalar aic_adj = -2*e(ll) + 2*e(k)

	estimates store T2M2
    
		/* This is the alternate version of this specification that we mention that omits sevvio. */
	
ml model lf wblreg (trgterra = logrelcap dyaddem3 contig rival ethnic terr socun) /ln_p, robust

	ml maximize

		/* Note that with no covariates the Newton-Raphson MLE with the startvalues option */
		/* produces a maximum at the largest observed duration. Graphing the hazards shows that */
		/* it has effectively flipped the two groups with the SS cdf reaching 1 around t=50 and */
		/* the W cdf reaching 1 near 1400. The others all find maxima at 44 or 20 depending on */
		/* the starting values (I manually entered them in exploring this). 44 is the best, though. */
	
warofatt trgterra, robust startvalues tstarmax(499) iterate(100) difficult

	estadd scalar aic_adj = -2*e(ll) + 2*(e(k) + 1)

	estimates store T2M3
	estimates save bds2020psrm-T2M3, replace
		
warofatt trgterra sevvio logrelcap dyaddem3 contig rival ethnic terr socun, robust startvalues tstarmax(499) ///
	iterate(100) difficult 

	estadd scalar aic_adj = -2*e(ll) + 2*(e(k) + 1)

	estimates store T2M4

		/* This is the alternate version of this specification that we mention that omits sevvio. */
	
warofatt trgterra logrelcap dyaddem3 contig rival ethnic terr socun, robust startvalues tstarmax(499) ///
	iterate(100) difficult 
		
		/* Include covariates in both equations. */
	
warofatt trgterra sevvio logrelcap dyaddem3 contig rival ethnic terr socun, probss(logrelcap dyaddem3) ///
	robust startvalues tstarmax(499) iterate(100) difficult 
		
	estadd scalar aic_adj = -2*e(ll) + 2*(e(k) + 1)

	estimates store T2M5
	estimates save bds2020psrm-T2M5, replace

		/* This is the alternate version of this specification that we mention that omits sevvio. */
	
warofatt trgterra logrelcap dyaddem3 contig rival ethnic terr socun, probss(logrelcap dyaddem3) ///
	robust startvalues tstarmax(499) iterate(100) difficult 
		
	
  	/**********************************************************/
	/* Do chi2 tests for model fit.	Note that the critical    */
	/* value needs to be adjusted via the chi-bar correction. */
  	/**********************************************************/


lrtest T2M3 T2M1, force
lrtest T2M4 T2M2, force
lrtest T2M5 T2M2, force
    
    
   	/*******************/
	/* Create table 2. */
   	/*******************/


estout T2M1 T2M2 T2M3 T2M4 T2M5 ///
	using bds2020psrm-table02.txt, ///
	cells(b(star fmt(3)) se(par)) ///
	starlevels(* 0.10 ** 0.05) ///
	equations(1:1:1:1:1 , 2:2:4:4:4) ///
	eqlabels("Weak Hazard" "Weak Dur. Dep." "SS Hazard" "SS Dur. Dep." "Probability Cured") ///
	order(#1: #2: StrongStrong: ln_p_SS: logit_pi:) ///
	stats(N ll tstarmax Fwtstar aic_adj, fmt(0 2 0 4 2 2)) ///
	label varwidth(30) modelwidth(12) ///
	varlabels(_cons constant) ///
	style(tab) replace

log close
clear
exit, STATA
